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Abstract: We compute the next-to-next-to-leading order (NNLO) QCD corrections to the 
six most important event shape variables related to three-particle final states in electron- 
positron annihilation. The corrections are sizeable for all variables, however their magni- 
tude is substantially different for different observables. We observe that the NNLO correc- 
tions yield a considerably better agreement between theory and experimental data both in 
shape and normalisation of the event shape distributions. The renormalisation scale de- 
pendence of the theoretical prediction is substantially reduced compared to the previously 
existing NLO results. Our results will allow a precise determination of the strong coupling 
constant from event shape data collected at LEP. 
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1. Introduction 

For more than a decade experiments at LEP (CERN) and SLC (SLAC) gathered a wealth 
of high precision high energy hadronic data from electron-positron annihilation at a range 
of centre-of-mass energies [1-5]. This data provides one of the cleanest ways of probing our 
quantitative understanding of QCD. This is particularly so because the strong interactions 
occur only in the final state and are not entangled with the parton density functions as- 
sociated with beams of hadrons. As the understanding of the strong interaction, and the 
capability of making more precise theoretical predictions, develops, more and more strin- 
gent comparisons of theory and experiment are possible, leading to improved measurements 
of fundamental quantities such as the strong coupling constant [6]. 

In addition to measuring multi-jet production rates, more specific information about 
the topology of the events can be extracted. To this end, many variables have been in- 
troduced which characterise the hadronic structure of an event. For example, we can ask 



- 1 - 



how planar or how collimated an event is. In general, a variable is described as n jet-like 
if it vanishes for a final state configuration of (n — 1) hadrons 1 . With the precision data 
from LEP and SLC, experimental distributions for such event shape variables have been 
extensively studied and have been compared with theoretical calculations based on next-to- 
leading order (NLO) parton-level event generator programs [7-9], improved by resumming 
kinematically-dominant leading and next-to-leading logarithms (NLO+NLL) [10, 11] and 
by the inclusion of non-perturbative models of power-suppressed hadronisation effects [12]. 

Comparing the different sources of error in the extraction of a s from hadronic data, 
one finds that the purely experimental error is negligible compared to the theoretical un- 
certainty. There are two sources of theoretical uncertainty: the theoretical description of 
the parton-to-hadron transition (hadronisation uncertainty) and the uncertainty stemming 
from the truncation of the perturbative series at a certain order, as estimated by scale 
variations (perturbative or scale uncertainty). Although the precise size of the hadroni- 
sation uncertainty is debatable and perhaps often underestimated, it is conventional to 
consider the scale uncertainty as the dominant source of theoretical error on the precise 
determination of a s from three-jet observables. 

For the bulk of the paper we are concerned with the next-to-next-to-leading order 
(NNLO) perturbative corrections to three jet-like shape variables. To be precise, we present 
the NNLO coefficients for the differential distributions of thrust, the wide and total jet 
broadening, heavy hemisphere mass, C parameter and the jet transition variable Y3 for 
the Durham jet algorithm. These results are obtained using a numerical implementation 
of the two-loop 7* — ► 3 parton [13], the one-loop 7* — ► 4 partons [14] and the tree-level 
7* — > 5 parton matrix elements [15]. Each of the contributions becomes singular when one 
or more partons are soft and/or collinear. In previous work, we have developed an antenna 
subtraction method [16] for isolating singularities and ensuring that the final result is 
infrared finite [17]. The resulting numerical program, EERAD3, yields the full kinematical 
information on the partonic final state and can be applied to generic infrared safe three-jet 
observables. 

In section 2, we provide definitions of the relevant three-jet shape variables while 
section 3 reviews the structure of the perturbative predictions. Section 4 gives a brief 
description of the NNLO calculation and its implementation in the multi-purpose parton 
level Monte Carlo program EERAD3. Results for the event shape distributions are reported 
in sections 5 and 6, together with an estimate of the remaining perturbative uncertainty due 
to variations of the renormalisation scale. The parton level predictions are also compared 
with hadron-level experimental data. Finally, our results are summarised in section 7. 

2. Event shape variables 

In order to characterise hadronic final states in electron-positron annihilation, a variety 
of event shape variables have been proposed in the literature, for a review see e.g. [18]. 

x It should be noted that sometimes in the literature, especially in works on resummation, event shapes 
requiring three particles are called two-jet event shapes, while those requiring four particles are called 
three-jet event shapes. 
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These variables can be categorised into two classes, according to the minimal number of 
final-state particles required for them to be non-vanishing: the most common variables 
require three particles (and are thus closely related to three-jet final states), while several 
other variables were constructed such that they require at least four particles (related to 
four-jet final states). 

Among the event shapes requiring three-particle final states, six variables were studied 
in great detail: the thrust T [19], the normalised heavy jet mass M 2 H js [20], the wide 
and total jet broadenings By/ and [21], the C-parameter [22] and the transition from 
three-jet to two-jet final states in the Durham jet algorithm Y3 [23]. 

(a) Thrust, T [19] 

The thrust variable for a hadronic final state in e + e~ annihilation is defined as [19] 

where pi denotes the three- momentum of particle i, with the sum running over all 
particles. The unit vector n is varied to find the thrust direction ut which maximises 
the expression in parentheses. 

The maximum value of thrust, T — > 1, is obtained in the limit where there are only 
two particles in the event. For a three-particle event the minimum value of thrust is 
T = 2/3. 

(b) Heavy hemisphere mass, M 2 H /s [20] 

In the original definition [20] one divides the event into two hemispheres. In each 
hemisphere, Hi, one also computes the hemisphere invariant mass as: 
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where E v [ s is the total energy visible in the event. In the original definition, the 
hemisphere is chosen such that M 2 + M| is minimised. We follow the more customary 
definition whereby the hemispheres are separated by the plane orthogonal to the 
thrust axis. 

The larger of the two hemisphere invariant masses yields the heavy jet mass: 

p = M 2 H /s = max(M 1 2 /s, M 2 2 /s) . (2.3) 

In the two-particle limit p — > 0, while for a three-particle event p < 1/3. 
The associated light hemisphere mass, 

Ml/a = min(M 1 2 /s, M 2 2 /s) (2.4) 

is an example of a four-jet observable and vanishes in the three-particle limit. 

At lowest order, the heavy jet mass and the (1 — T) distribution are identical. How- 
ever, this degeneracy is lifted at next-to-leading order. 
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(c) Jet Broadening, By/ and [21] 

Taking a plane perpendicular to ny through the coordinate origin, one defines two 
event hemispheres B\^- In each of them, one determines the hemisphere broadening: 



^2 \fk x n T \ 

Mi 



Bi = keH * . (2.5) 



k 

The wide and total jet broadening are then defined as 

B w = max(B 1 ,B 2 ) , (2.6) 
B T = Bi + B 2 . (2.7) 

In the two-particle limit By/ — > and St — > 0. The maximum broadening for a 
three-particle event is 5t = -Bvy = l/(2\/3). 

The narrow jet broadening, 

B N = mm(B 1 ,B 2 ) , (2.8) 

is another four-jet observable and vanishes when only three particles are in the event. 

(d) The C parameter, [22] 

The linearised momentum tensor 

Q^ = —1—\2t^T' (a, = 1,2,3), (2.9) 

has three eigenvalues Aj, which are used to construct the C-parameter: 

C = 3 (AiA 2 + A 2 A 3 + A 3 Ai) . (2.10) 
This definition is equivalent to 

c = 3 (e n e 22 + e 22 e 33 + e 33 e n - e 12 e 12 - e 23 e 23 - e 31 e 31 ) . (2.11) 

The related four-jet observable is the L>-parameter, 

£> = 27AiA 2 A 3 . (2.12) 

(e) The jet transition variable, Y3 [23] 

The jet transition variable I3 is defined as the value of the jet resolution parameter 
y cu t for which an event changes from a three-jet to a two-jet configuration with some 
jet defining scheme. 
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Here, we focus on the Durham jet algorithm which clusters particles into jets by 
computing the measurement variable 



2min(£ t 2 ,£|)(l-cos%) 
VH,D = y 2 (2-13) 

vis 

for each pair of particles. The pair with the lowest Uij,D is replaced by a 

pseudoparticle whose four-momentum is given by the sum of the four-momenta of 
particles i and j ('E' recombination scheme). This procedure is repeated as long as 
pairs with invariant mass below the predefined resolution parameter yij t D < Vcut are 
found. Once the clustering is terminated, the remaining (pseudo-)particles are the 
jets. 



3. Event shapes in perturbation theory 

The perturbative expansion for the distribution of a generic observable y up to NNLO at 
centre-of-mass energy y/s for renormalisation scale fi 2 = s and a s = a s (y / i) is given by 

1 do- fa s \ dA fa s \2 dB /a s \3 dC 4 

Here the event shape distribution is normalised to the total hadronic cross section Chad- 
With the assumption of massless quarks, then at NNLO we have, 

= •» i 1 + \ c - (£) + (S) 2 + ) • < 3 ' 2 ' 

where the Born cross section for e + e~ — ► qq is 

ao = — iVe*. (3.3) 



The constant is given by, 



3 / 123 

-^fl + C F C A ( — - 44C 3 ) + C F T R N F (-22 + I6C3) 



(3.4) 



where the QCD colour factors are 



N 2 — I 1 
C A = N, C F = ^ r , T R = - (3.5) 

for N = 3 colours and 7Vp light quark flavours. 

In practice, we compute the perturbative coefficients A, B and C, which are all nor- 
malised to uq: 

1 da /a.,\ dA /a. s \ 2 d-B /a s \ 3 dC 4 
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However, A, B and C are straightforwardly related to A, B and C, 

A = A , 
B = B 



\C F A, 



C 



C--C F B + 



* n 2 

7°f 



(3.7) 



These coefficients are computed at a renormalisation scale fixed to the centre-of-mass en- 
ergy, and depend therefore only on the value of the observable y. They explicitly include 
only QCD corrections with non-singlet quark couplings and are therefore independent of 
electroweak couplings. At 0(a 2 ), these amount to the full corrections, while the 0(a^) 
corrections also receive a pure-singlet contribution. This pure-singlet contribution arises 
from the interference of diagrams where the external gauge boson couples to different quark 
lines. In four-jet observables at O(af), these singlet contributions were found to be ex- 
tremely small [24]. Also, the pure-singlet contribution from three-gluon final states to 
three-jet observables was found to be negligible [25]. This small correction to NNLO is 
denoted by 5c- 



1 da 

Co dy 



d5 c 
dy 



(s, Mz, a, sin 2 Oyy, c q 



(3.8) 



NNLO, pure singlet 

where c q denotes the set of all electroweak vector and axial-vector quark couplings. 

First-order electroweak corrections to event shape observables could be of a magnitude 
comparable to the NNLO QCD corrections. Like the pure-singlet NNLO contributions, 
these do also not factorise onto Co- The first-order electroweak corrections affect the dis- 
tribution itself and the normalisation <7had- Collectively, they give a contribution of the 
form, 



1 do- 



c-had dy 



/ a 



\2irJ V2vr, 



d5E\v 
dy 



(s,M z , a, sin 2 @ w ,c q ). 



(3.9) 



electroweak, 0(aa s ) 

These corrections are not complete at present [26], and clearly deserve further study. 

In summary, the expression for event shape distributions accurate to NNLO in QCD 
and NLO in the electroweak theory reads: 



1 do- 



— - f^l) — + f^lY — + ( as Y dU 

dv ~ V2tt/ dy V2vr/ dv \2ttJ dv 



3 dC 



Ohad dy 



,2vr 
/a 



dy \2irJ dy \2irJ dy 



\ 3 d5c o 
_) —(s,M z ,a,sm^e w ,c q ) 



(3.10) 



In the following, we will focus on the QCD non-singlet expression (3.1), since 5c can be 
safely neglected, and the computation of 5ew needs further work. 

The QCD coupling constant evolves according to the renormalisation group equation, 
which is to NNLO: 



> da s (n) 
d/j, 2 



-a s (n) 



A) 



« s (/j) 
2vr 



Q s (/i) 

2vr 



+ P2 



a s (p) 
2vr 



+ 0(a 



(3.11) 
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with the MS-scheme coefficients 
11C a -4TrN f 



Po = 
Pi 



6 

11C\ - 10CaT r N f - QC F T R N F 



6 

02 = ^(2857^1 + W8C 2 F T R N F - 1230C F C A T R N F - 2830C 2 A T R N F 

+264C F T^N F + 316C A T R N F ) . (3.12) 

Equation (3.11) is solved by introducing A as integration constant with L = log( / u 2 /A 2 ), 
yielding the running coupling constant: 



27T / ftlQgL 1 (0l ( 2 l^-M^ /q i q\ 

a * (M) = Hz. v - $$ — + w? U (log L " logL " 1} + ) ■ (3 - 13) 



In terms of the running coupling a s (fi), the NNLO (non-singlet) expression for event 
shape distributions becomes 

1 da. 2 . fa 9 (n)\dA fa s (fx)\ 2 fdB dA a , ^ 2 



a had dy \ 2tt J dy \ 2vr / \ dy dy s 



+ -r^ hr + 2 T- A) log — + T" /?o log - + /?ilog 

Z7r / \ dy dy s dy \ s s 



+C?(<4) • (3.14) 



4. Calculation of NNLO corrections 



Three-jet production at tree-level is induced by the decay of a virtual photon (or other 
neutral gauge boson) into a quark-antiquark-gluon final state. At higher orders, this process 
receives corrections from extra real or virtual particles. The individual partonic channels 
that contribute through to NNLO are shown in Table 1. All of the tree- level and loop 
amplitudes associated with these channels are known in the literature [13-15,27]. 

For a given partonic final state, the event shape observable y is computed according to 
the same definition as in the experiment, which is applied to partons instead of hadrons. 
At leading order, all three final state partons must be well separated from each other, such 
that y differs from the trivial two-parton limit. At NLO, up to four partons can be present 
in the final state, two of which can be clustered together, whereas at NNLO, the final state 
can consist of up to five partons, and as many as three partons can be clustered together. 
The more partons in the final state, the better one expects the matching between theory 
and experiment to be [28]. 

The two-loop 7* — ► qqg matrix elements were derived in [13] by reducing all relevant 
Feynman integrals to a small set of master integrals using integration-by-parts [29] and 
Lorentz invariance [30] identities, solved with the Laporta algorithm [31]. The master 
integrals [32] were computed from their differential equations [30] and expressed analytically 
in terms of one- and two-dimensional harmonic polylogarithms [33]. 
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LO 7* — > qqg tree level 

NLO 7* — ► g qg one loop 

7* ~^ QQ99 tree level 

7* —> qqqq tree level 

NNLO 7* — >• qqg two loop 

7* — ► qqgg one loop 

7*^qqqq one loop 

7* ~^ QQQQ9 tree level 

1* ^ 1Q9 9 9 tree level 

Table 1: Non-singlet partonic contributions to three-jet event shape observables in perturbative 
QCD. 



The one-loop four-parton matrix elements relevant here [14] were originally derived in 
the context of NLO corrections to four-jet production and related event shapes [34,35]. 
One of these four-jet parton-level event generator programs [35] is the starting point for 
our calculation, since it already contains all relevant four-parton and five-parton matrix 
elements. 

The four-parton and five-parton contributions to three-jet-like final states at NNLO 
contain infrared real radiation singularities, which have to be extracted and combined 
with the infrared singularities [36] present in the virtual three-parton and four-parton 
contributions to yield a finite result. In our case, this is accomplished by introducing 
subtraction functions, which account for the infrared real radiation singularities, and are 
sufficiently simple to be integrated analytically. Schematically, this subtraction reads: 



da nnlo = / (dcrNNLO ~ da NNLO ) 



d$ 5 

V,l 

NNLO KXU NNLO 



+ I [daX}\iT.n ~ da Vb ' 

d<E>4 



V,2 

NNLO ' 



+ / da NNLO + / da^ S / LO + / da) 

Jd&c, Jd&4 Jd&s 

where da NNLO denotes the real radiation subtraction term coinciding with the five-parton 
tree level cross section da^ NLO in all singular limits [37]. Likewise, da^ NLO is the one- loop 
virtual subtraction term coinciding with the one-loop four-parton cross section da^ l NLO in 
all singular limits [38] . Finally, the two-loop correction to the three-parton cross section is 
denoted by da^ NL0 . With these, each line in the above equation is individually infrared 
finite, and can be integrated numerically. 

Systematic methods to derive and integrate subtraction terms were available in the 
literature only to NLO [39,40]. Physical results for the special case of NNLO Higgs pro- 
duction have been achieved in [41]. In the context of this project, we fully developed an 
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NNLO subtraction formalism [16,17,42], based on the antenna subtraction method orig- 
inally proposed at NLO [35,40]. The basic idea of the antenna subtraction approach is 
to construct the subtraction terms from antenna functions. Each antenna function encap- 
sulates all singular limits due to the emission of one or two unresolved partons between 
two colour-connected hard partons. This construction exploits the universal factorisa- 
tion of phase space and squared matrix elements in all unresolved limits. The individual 
antenna functions are obtained by normalising three-parton and four-parton tree-level ma- 
trix elements and three-parton one-loop matrix elements to the corresponding two-parton 
tree-level matrix elements. Three different types of antenna functions are required, cor- 
responding to the different pairs of hard partons forming the antenna: quark-antiquark, 
quark-gluon and gluon-gluon antenna functions. All these can be derived systematically 
from matrix elements [42] for physical processes. 

The factorisation of the final state phase space into antenna phase space and hard phase 
space requires a mapping of the antenna momenta onto reduced hard momenta. We use the 
mapping derived in [43] for the three-parton and four-parton antenna functions. To extract 
the infrared poles of the subtraction terms, the antenna functions must be integrated 
analytically over the appropriate antenna phase spaces, which is done by reduction [44] to 
known phase space master integrals [45]. 

We tested the proper implementation of the subtraction by generating trajectories of 
phase space points approaching a given single or double unresolved limit. Along these 
trajectories, we observe that the antenna subtraction terms converge towards the physical 
matrix elements, and that the cancellations among individual contributions to the subtrac- 
tion terms take place as expected. Moreover, we checked the correctness of the subtraction 
by introducing a lower cut (slicing parameter) on the phase space variables, and observing 
that our results are independent of this cut (provided it is chosen small enough). This 
behaviour indicates that the subtraction terms ensure that the contribution of potentially 
singular regions of the final state phase space does not contribute to the numerical inte- 
grals, but is accounted for analytically. A detailed description of the calculation can be 
found in [17]. 

The resulting numerical program, EERAD3, yields the full kinematical information on a 
given multi-parton final state. It can thus be used to compute any infrared-safe observable 
related to three-particle final states at 0(o%.) in e + e~-annihilation. 

5. NNLO distributions 

In this section, we discuss the size and shape of the LO, NLO and NNLO coefficients of 
the perturbative expansion of the various event shape observables defined in Eq. (3.6). For 
convenience, we weight the distribution by the observable. 

The precise size and shape of the NNLO corrections depend on the observable in 
question. However, all contributions are dominated by the behaviour in the two-jet region 
where the observable generally tends to zero. Of course, typical hadronic events contain 
many hadrons and it is extremely unlikely that the value of any event shape is precisely 
zero for any experimental event. However, in the fixed order partonic calculation, where 
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there are at most five partons present in the final state, one or more of the partons may 
be soft and/or collinear, and the observable may approach zero. In such circumstances, 
soft gluon singularities cause the fixed order prediction to become wildly unstable and grow 
logarithmically. In the infrared limit y — > 0, the perturbative coefficients have the following 
form, 

d ^ ATA 

y- AxL + Aq 

Ay 

y B 3 L 3 + B 2 L 2 + B\L + B 

ay 

y^- ~ C 5 L 5 + C 4 L 4 + C 3 L 3 + C 2 L 2 + C ± L + C (5.1) 
dy 

where L = ln(l/y) and C n are (as yet) undetermined coefficients. Whenever L is sufficiently 
large, resummation effects will be important. In our numerical studies, we therefore impose 
a cut on the size of y which is typically in the range 0.001 - 0.01, since for such small values 
of y we do not trust the fixed order prediction. 

Even away from the infrared region, the shape of the fixed order prediction is heavily 
influenced by cancellations between the real and virtual contributions. The LO contribution 
A is very large and positive at small y and decreases monotonically as y increases. The 
NLO contribution B is negative at small y, but exhibits a turn-over, typically at y ~ 0.05. 
Similarly, the NNLO contribution C also exhibits a turn-over, but at a slightly larger value 
of y. The precise positions of the maxima of the distributions depend on the observable 
under consideration. 

A second generic feature occurs when the paucity of final state particles imposes a 
maximum value for the observable. Examples include (1 — T) and C which are required 
to be less than 0.33 and 0.75 respectively for three-parton final states. As the number of 
partons increases with the perturbative order, this limit is relaxed and larger values of the 
observable are accessed. 

Finally, typical values of the strong coupling constant lie around a s ~ 0.12, so that 
^ ~ 1/50. It is well known that the NLO corrections are large, B y ~ (15 — 30)^4^, leading 
to a 30-60% NLO effect in the region where the perturbative calculation is expected to be 
reliable. However, we observe that in all cases, the NNLO coefficients are also significant, 
C y ~ (200 - 800) A y , leading to a further 7-28% NNLO correction. 

5.1 Thrust 

Thrust is defined in section 2(a). First results for the NNLO corrections to the thrust 
distribution were presented in Ref. [46]. The perturbative coefficients for the thrust dis- 
tribution weighted by (1 — T) are shown in Fig. 1. As discussed earlier, the shape of 
the contribution is dominated by the infrared region at (1 — T) — ► 0. At small (1 — T), 
the LO contribution A is very large and positive, while the NLO and NNLO coefficients 
B and C are rising and exhibit a turn-over at moderate values of (1 — T). We observe 
that the peak moves from about 0.04 (NLO) to 0.06 (NNLO). We also see that the NLO 
and NNLO distributions progressively extend to larger and larger values of (1 — T) as the 
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Figure 1: Coefficients of the leading order, next-to- leading order and next-to-next-to- leading order 
contributions to the thrust distribution as defined in Eq. (3.6) and weighted by (1 — T). The dotted 
line in the C coefficient indicates the distribution prior to correction of the soft large-angle radiation 
terms (see erratum at the end of the paper). 




0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 



P P P 

Figure 2: Coefficients of the leading order, next-to-leading order and next-to-next-to-leading order 
contributions to the heavy jet mass distribution as defined in Eq. (3.6) and weighted by p. The 
dotted line in the C coefficient indicates the distribution prior to correction of the soft large-angle 
radiation terms. 

phase space restrictions on large values of (1 — T) are relaxed. In the intermediate region, 
0.04 < (1 — T) < 0.33, we observe that the perturbative coefficients are roughly in the 
ratio, A : B : C ~ 1 : 30 : 800. Setting a s ~ 0.12 and using Eq. (3.7), this indicates 
corrections which are of relative magnitude LO : NLO : NNLO ~ 1 : 0.53 : 0.27, such that 
the NNLO corrections increase the NLO prediction by another 18%. 

5.2 Heavy jet mass 

The definition of the heavy jet mass given in section 2(b) is the larger invariant mass of 
the two hemispheres formed by separating the event by a plane normal to the thrust axis. 
The perturbative coefficients for the heavy jet mass distribution weighted by p are shown 
in Fig. 2. At lowest order, the heavy jet mass and the (1 — T) distribution are identical, 
so that A does not extend past p = 0.33. At higher orders, the distribution extends to 
larger values, with a small negative NNLO contribution around 0.33. In the intermediate 
region, 0.02 < p < 0.33, the perturbative coefficients are roughly A : B : C ~ 1 : 20 : 400 
indicating corrections of approximately LO : NLO : NNLO ~ 1 : 0.34 : 0.13, translating into 
a 10% enhancement of NNLO over NLO. Comparing Fig. 1(b) with 2(b) and Fig. 1(c) with 
2(c) we see clearly the rather different behaviour of the higher order corrections to these 
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0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 



Bj 

Figure 3: Coefficients of the leading order, next-to-leading order and next-to-next-to-leading order 
contributions to the total jet broadening distribution as defined in Eq. (3.6) and weighted by Bt- 
The dotted line in the C coefficient indicates the distribution prior to correction of the soft large- 
angle radiation terms. 




Figure 4: Coefficients of the leading order, next-to-leading order and next-to-next-to-leading order 
contributions to the wide jet broadening distribution as defined in Eq. (3.6) and weighted by By/. 
The dotted line in the C coefficient indicates the distribution prior to correction of the soft large- 
angle radiation terms. 

observables, particularly in the region beyond the LO kinematic bound where partonic 
configurations with two or more partons in each hemisphere contribute differently to each 
observable. 

5.3 Jet broadenings 

The jet broadenings are defined in section 2(c) by dividing the event into two hemispheres 
using a plane normal to the thrust axis.. At lowest order By/ and Bt are identical, but 
their distributions receive different higher order corrections from partonic configurations 
with two or more partons in each hemisphere. 

The perturbative coefficients for the Bt (Bw) distributions weighted by Bt {Bw) 
are shown in Fig. 3 (Fig. 4) respectively. The structures evident around Bt, Bw ~ 
(l/2\/3) ~ 0.29 are generated by four and five parton events and are therefore different for 
the two observables. For more moderate Bt values between 0.04 and 0.29, the perturbative 
coefficients are in the ratio A : B : C ~ 1 : 35 : 800. Including the factors of a s , this 
leads to corrections LO : NLO : NNLO ~ 1 : 0.63 : 0.27 for a s ~ 0.12, which amounts 
to NNLO corrections of 17% of the NLO result. We observe that the corrections for By/ 
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Figure 5: Coefficients of the leading order, next-to-leading order and next-to-next-to-leading order 
contributions to the C parameter distribution as defined in Eq. (3.6) and weighted by C. The dotted 
line in the C coefficient indicates the distribution prior to correction of the soft large-angle radiation 
terms. 

are considerably smaller than those for Bt, A : B : C ~ 1 : 20 : 400 or equivalently 
LO : NLO : NNLO ~ 1 : 0.34 : 0.13, which yields a 10% NNLO effect over NLO. When 
Bt < 0.04 (B\y < 0.04), infrared logarithms must be resummed to produce a meaningful 
result. 

5.4 C-parameter 

The C parameter is defined in section 2(d) and the perturbative distributions at LO, NLO 
and NNLO weighted by C are shown in Fig. 5. The LO kinematic limit at C = 0.75 is 
clearly visible. At NLO (and NNLO), four (and five) parton events can generate larger 
values of C, leading to a sharp peak around C ~ 0.75. The approximate size of the 
corrections for 0.1 < C < 0.75 is A : B : C ~ 1 : 30 : 700, or, including the factors 
of (a s /27r) with a s ~ 0.12, in the ratio LO : NLO : NNLO ~ 1 : 0.53 : 0.23, resulting 
in a 15% enhancement of NNLO over NLO. At smaller values of C < 0.1, large infrared 
logarithms render the fixed order prediction unreliable and must be resummed. Similarly, 
large logarithms are produced around the LO kinematic limit, C ~ 0.75 which must also 
be resummed. 

5.5 y 3 

The jet transition variable I3 is defined in section 2(e). It describes the value of the 
jet resolution parameter y cut for which an event changes from a three-jet to a two-jet 
configuration within the Durham jet algorithm. The perturbative distributions at LO, NLO 
and NNLO weighted by I3 are shown in Fig. 5. As with all of the event shapes, Y^dA/dY^ 
is linear when plotted on a logarithmic scale. For moderate values of I3, 2 < — ln(Y3) < 6, 
the corrections are positive. In this region, the approximate size of the corrections is 
A : B : C ~ 1 : 15 : 200, or, including the factors of (a s /27r) with a s ~ 0.12, in the ratio 
LO : NLO : NNLO ~ 1 : 0.25 : 0.06, which produces a 5% NNLO effect over NLO. However, 
at smaller values of I3 (larger values of -ln(l3)) resummation of logarithmic contributions 
are clearly mandatory. 
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Figure 6: Coefficients of the leading order, next-to-leading order and next-to-next-to-leading order 
contributions to the distribution of the jet transition variable Y3 as defined in Eq. (3.6) and weighted 
by Y 3 . The dotted line in the C coefficient indicates the distribution prior to correction of the soft 
large-angle radiation terms. 



6. Comparison with data 

We have presented the NNLO corrections to six event-shape distributions. As we have 
shown in the previous section, the magnitude of the NNLO correction is different for the 
six variables. 

Each of the event shapes considered here has been studied in depth by all four ex- 
periments at LEP at centre-of-mass energies of 91.2, 133, 161, 172, 183, 189, 200 and 206 
GeV [1-4]. Within the experimental uncertainties, these data sets are mutually consistent. 

In this paper, we select data from ALEPH [1] as a representative set of hadronic final 
states in electron-positron annihilation to illustrate the improvement in the theoretical 
prediction due to the inclusion of the NNLO perturbative contribution. The only free 
parameter in our predictions is the strong coupling constant; we use the current world 
average value a s {M z ) = 0.1189 [47]. 

The experimental event-shape distributions were computed using the reconstructed 
momenta and energies of charged and neutral particles. The measurements have been cor- 
rected for detector effects and the final distributions correspond to the particle (or hadron) 
level (stable hadrons and leptons after hadronisation). In addition, at LEP2 energies above 
the Z peak the data were corrected for initial-state radiation effects and backgrounds, 
mainly from VF-pair production, were subtracted. The experimental uncertainties were 
estimated by varying event and particle selection cuts and are below 1% at LEP1 and 
between 0.5% and 1.5% at LEP2. For further details we refer the interested reader to 
Ref. [1]. 

6.1 Thrust 

Figure 7 displays the perturbative expression for the thrust distribution 2 at LO, NLO and 
NNLO, evaluated at Q = Mz- The error band indicates the variation of the prediction 
under shifts of the renormalisation scale in the range £ [Q/2;2Q] around the e + e~ 
centre-of-mass energy Q. The relative scale uncertainty is reduced by about 30% between 
NLO and NNLO. 

2 First results for the NNLO corrections to the thrust distribution were presented in Ref. [46] 
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Figure 7: Thrust distribution at Q = M z at LO (blue), NLO (green) and NNLO (red). The solid 
lines represent the prediction for renormalisation scale /i = Q and a s (Mz) = 0.1189, while the 
shaded region shows the variation due to varying the renormalisation scale between [i = Q/2 and 
(X = 2Q. The data is taken from [1]. 



The inclusion of the NNLO corrections enhances the thrust distribution by around 
(15-20)% over the range 0.04 < (1 — T) < 0.33, where — ln(l — T) is not too large. Outside 
this range, one does not expect the perturbative fixed-order prediction to yield reliable 
results. For (1 — T) — * 0, the convergence of the perturbative series is spoilt by powers 
of logarithms ln(l — T) appearing in higher perturbative orders, thus necessitating an all- 
order resummation of these logarithmic terms [10, 11], and a matching of fixed-order and 
resummed predictions [48]. 

The perturbative parton-level prediction is compared with the hadron-level data from 
the ALEPH collaboration [1] in Figure 7 and Figure 8. We observe that for all Q values, 
the shape and normalisation of the parton level NNLO prediction agrees better with the 
data than at NLO. We also see that the NNLO corrections account for approximately half 
of the difference between the parton-level NLO prediction and the hadron-level data. 

6.2 Heavy jet mass 

The perturbative prediction for the heavy jet mass distribution is displayed in Figure 9. 
The solid lines represent the prediction at the physical scale Q = Mz, while the shaded 
bands represent the effect of varying the renormalisation scale upwards and downwards 
by a factor of 2. We observe that the relative scale uncertainty is reduced by about 50% 
between NLO and NNLO. It is noteworthy that the original motivation for introducing the 
heavy jet mass distribution [20] was the hope for improved perturbative stability over the 
thrust distribution. This improved stability was not evident from the existing NLO results 
alone, but becomes visible at NNLO. 

Compared to NLO, the inclusion of the NNLO corrections enhances the heavy jet 
mass distribution by around 10% over the range 0.02 < p < 0.33, where ln(p) is not too 
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Figure 8: The thrust distribution (with fj, = Q and a s (M z ) = 0.1189) at LO (dotted), 
NLO (dashed) and NNLO (solid) compared to experimental data from ALEPH [1] for Q = 
133 GeV, . . . , 206 GeV. 



- 16 - 



"D 



0.6 



Q_ 



T3 
CO 



0.1 



0.4 



0.3 



0.5 



0.2 








a s (M z ) = 0.1189 







0.1 



0.2 



0.3 



0.4 



P 



Figure 9: Heavy jet mass distribution at Q = Mz at LO (blue), NLO (green) and NNLO (red). 
The solid lines represent the prediction for renormalisation scale \i = Q and a s (Mz) — 0.1189, while 
the shaded region shows the variation due to varying the renormalisation scale between p = Q/2 
and p = 2Q. The data is taken from [1]. 

large. At smaller p values, large ln(l/p) corrections must be resummed to all orders [49] 
and matched onto the perturbative prediction. Nevertheless, in the moderate to large p 
region, the NNLO corrections render the fixed order prediction significantly closer to the 
experimental data [1]. 

Figure 10 shows the prediction for a range of Q values together with the hadron-level 
data from the ALEPH collaboration [1]. For this observable, the NNLO corrections are 
relatively small, however, for all Q values, the shape and normalisation of the parton-level 
NNLO prediction agrees slightly better with the hadron-level data than at NLO. 

6.3 Jet broadenings 

Predictions for the total and wide jet broadenings are displayed in Figures 11 and 12. 
The solid lines represent the prediction at the physical scale Q = Mz, while the shaded 
bands represent the effect of varying the renormalisation scale upwards and downwards by 
a factor of 2. We observe that the relative scale uncertainty in the Bt (By/) distribution 
is reduced by about 40% (50%) between NLO and NNLO. 

As anticipated from the discussion in section 5.3, we observe that the perturbative 
corrections are smaller for By/ than for Bt- In the region where perturbation theory 
is expected to yield reliable results, (Bt, By/) > 0.05, we observe an enhancement of 
(15-20)% in Bt and of (8-12)% in By/. As with (1 — T) and the heavy jet mass, the 
two broadenings are identical at leading order, but display a largely different behaviour 
in the higher perturbative corrections. At smaller values of broadening, large logarithmic 
corrections occur which must be resummed [21]. 
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Figure 10: Heavy jet mass distribution (with fi = Q and a s (Mz) = 0.1189) at LO (dot- 
ted), NLO (dashed) and NNLO (solid) compared to experimental data from ALEPH [1] for 
Q = 133 GeV, . . . , 206 GeV. 
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Figure 11: Total jet broadening distribution at Q — Mz at LO (blue), NLO (green) and NNLO 
(red). The solid lines represent the prediction for renormalisation scale /x = Q and a s (Mz) = 0.1189, 
while the shaded region shows the variation due to varying the renormalisation scale between 
/j = Q/2 and fj, = 2Q. The data is taken from [1]. 




Figure 12: Wide jet broadening distribution at Q = Mz at LO (blue), NLO (green) and NNLO 
(red). The solid lines represent the prediction for renormalisation scale [i — Q and a s (Mz) — 0.1189, 
while the shaded region shows the variation due to varying the renormalisation scale between 
/j = Q/2 and /i = 2Q. The data is taken from [1]. 
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Figure 13: Total jet broadening distribution (with fj, = Q and a s (Mz) = 0.1189) at LO 
(dotted), NLO (dashed) and NNLO (solid) compared to experimental data from ALEPH [1] for 
Q = 133 GeV, . . . , 206 GeV. 
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Figure 14: Wide jet broadening distribution (with fj, = Q and a s (Mz) = 0.1189) at LO 
(dotted), NLO (dashed) and NNLO (solid) compared to experimental data from ALEPH [1] for 
Q = 133 GeV, . . . , 206 GeV. 
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Figure 15: C parameter distribution at Q = Mz at LO (blue), NLO (green) and NNLO (red). The 
solid lines represent the prediction for renormalisation scale fj, = Q and a s (Mz) = 0.1189, while 
the shaded region shows the variation due to varying the renormalisation scale between [i = Q/2 
and fi = 2Q. The data is taken from [1]. 

To guide the eye, Figures 11 and 12 also show hadron-level data from the ALEPH 
collaboration [1]. For both broadenings, we see that the NNLO prediction lies closer to the 
data, and, in fact, accounts for much of the difference between the NLO prediction and the 
hadron-level data. 

The experiments at LEP also gathered data at higher Q values; Figures 13 and 14 com- 
pare the parton-level prediction at Q = 133 GeV, . . . , 206 GeV with hadron-level data from 
the ALEPH collaboration [1]. We observe that for all Q values, shape and normalisation 
of the parton level NNLO prediction agrees better with the data than at NLO. 

6.4 C-parameter 

The C parameter is one of the classic event shape observables and we display the perturba- 
tive prediction in Figure 15. The solid lines represent the prediction at the physical scale 
Q = Mz, while the shaded bands represent the effect of varying the renormalisation scale 
upwards and downwards by a factor of 2. We observe that the relative scale uncertainty 
is reduced by about 40% between NLO and NNLO. The NNLO corrections enhance the C 
parameter distribution by around (12-20)% over the range 0.1 < C < 0.75, where ln(l/C) 
is not too large. Figure 15 also shows hadron-level data from the ALEPH collaboration [1] 
and we observe that the NNLO parton-level prediction lies significantly closer to the data, 
and in fact, accounts for about one third of the difference between the NLO prediction and 
the data. 

At small C, one expects large logarithmic contributions ln(l/C) appearing in higher 
perturbative orders, thus necessitating an all-orders resummation of these logarithmic 
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Figure 16: C parameter distribution (with /i = Q and a s (Mz) = 0.1189) at LO (dot- 
ted), NLO (dashed) and NNLO (solid) compared to experimental data from ALEPH [1] for 
Q = 133 GeV, . . . , 206 GeV. 
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terms [50]. There are also large logarithms around C ~ 0.75, due to soft gluon divergences 
within the physical region (producing a so-called Sudakov shoulder in the distribution) 
which must also be resummed [51] to all orders. 

Figure 16 shows the prediction for a range of Q values together with the hadron-level 
data from the ALEPH collaboration [1]. For all Q values, the shape and normalisation of 
the parton level NNLO prediction agrees slightly better with the data than at NLO. 




2 4 6 8 10 



-ln(Y 3 ) 

Figure 17: The distribution for the jet transition variable, I3 at Q = Mz at LO (blue), NLO 
(green) and NNLO (red). The solid lines represent the prediction for renormalisation scale /i = Q 
and a s (Mz) = 0.1189, while the shaded region shows the variation due to varying the renormali- 
sation scale between /i = Q/2 and /i = 2Q. The data is taken from [1]. 

6.5 y 3 

Figure 17 displays the perturbative expression for the I3 distribution at LO, NLO and 
NNLO, evaluated at Q = Mz- The error band indicates the variation of the prediction 
under shifts of the renormalisation scale in the range [i E [Q/2;2Q] around the e + e~ 
centre-of-mass energy Q. The relative scale uncertainty is reduced by about 50% between 
NLO and NNLO. 

The NLO and NNLO corrections change the shape of the distribution considerably and 
introduce a turnover at — ln(Y3) ~ 5 — 6. We observe that the NNLO corrections modify 
the I3 distribution by around (3-5)% over the range 2 < — ln(Y~3) < 6, where — ln(l3) 
is moderate. At larger — ln(l3), one does not expect fixed-order perturbation theory to 
yield reliable results and the large infrared logarithms of the type a^ln^Y-j) must be 
resummed [52]. 

The perturbative parton-level prediction is compared with the hadron-level data from 
the ALEPH collaboration [1] in Figure 17 and Figure 18. We see that the quality of the 
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Figure 18: The Y 3 distribution (with fj, = Q and a s (M z ) = 0.1189) at LO (dotted), NLO (dashed) 
and NNLO (solid) compared to experimental data from ALEPH [1] for Q = 133 GeV, . . . , 206 GeV. 
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agreement between fixed order perturbation theory and data is much more Q dependent 
than for the other observables considered in this paper. At Q = Mz, the data is much 
more sharply peaked than the NNLO prediction. However, at higher Q values shown in 
Figure 18, the agreement between the NNLO prediction and the data around the peak 
region is very good - and significantly better than at NLO. 

7. Conclusions and Outlook 

The main goal of this paper is to provide improved theoretical predictions for hadronic 
event shape observables in electron-positron annihilation. To this end, we have presented 
new results for the next-to-next-to-leading order contributions to a number of important 
three-jet-like event shape observables in e + e~ collisions. These results are obtained using 
a numerical program, that is based on the matrix elements for 7* — > 3 partons at two-loop, 
7* — ► 4 partons at one-loop and 7* — ► 5 partons at tree-level. Each of these contributions 
becomes singular when one or more partons are soft and/or collinear, and we have developed 
and implemented an NNLO subtraction formalism to subtract these singularities, thereby 
yielding a finite NNLO prediction. The resulting numerical program, EERAD3, yields the full 
kinematical information on the partonic final state and can be applied to generic infrared 
safe three-jet observables. 

For the six event shapes considered here, and in kinematical regions where infrared 
logarithms are small enough to render their resummation unnecessary, the NLO corrections 
are generally large - of the order 30-60%. The NNLO effects produce a further 5-20% 
correction. Comparisons with existing data from LEP indicate an improved agreement with 
hadronic data and the fixed order NNLO parton-level prediction. In addition, the remaining 
theoretical uncertainty estimated by varying the renormalisation scale by a factor of two 
around the physical scale is also significantly reduced, typically by 30-50%. Importantly, 
the size of the corrections is different for different observables. 

Our results for the NNLO corrections open up a whole new range of possible compar- 
isons with the LEP data. For meaningful comparisons, one has to account for hadronisation 
effects, either by introducing hadron-level to parton-level correction factors, or by includ- 
ing power-suppressed hadronisation effects in the theoretical description. A first direct 
determination of the strong coupling constant from a fit of next-to-next-to-leading order 
QCD predictions to event-shape variables over a range of Q values will be reported in a 
separate publication [53] and should yield a much more precise value of a s (Mz) than that 
previously extracted from event shapes. Our predictions can be further improved by a NLL 
resummation of the large infrared logarithms that are present as y — > 0. The ingredients 
for lni? matching to NNLO are available in [11]. Studies in this direction are in progress 
and should yield a further improvement on the measurement of a 8 (Mz)- Similarly, our 
results will also allow a renewed study of power corrections, now matched to NNLO. 
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Erratum added 

Our implementation of NNLO corrections to three-jet-like observables [17], which was used 
for the present work, was checked by two subsequent studies: the calculation of all loga- 
rithmically enhanced contributions to the thrust distribution by Becher and Schwartz [A], 
and an independent implementation of our subtraction formulae by Weinzierl [B]. 

These works uncovered numerical discrepancies in the two-jet limit of the observables 
in two of the six colour factor contributions: N 2 and N°. In [B], it was shown that the 
origin of these discrepancies is in an oversubtraction of large-angle soft radiation. We 
described the corrected treatment of these terms in the erratum to [17]. 

As a consequence, the numercial values of the NNLO coefficients of all event shape 
distributions were modified. The new coefficients are displayed by the solid lines in Figures 
1-6, our original results are displayed there for comparison as dotted lines. It can be 
seen that in the genuine three-jet region, which is relevant for precision phenomenology, 
the changes have a minor numerical impact. The modification to the full event shape 
distributions is too small to be visible, except for I3 in the deep two-jet region. We 
therefore refrain from presenting revised figures 7-18. 
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